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NOMENCLATURE 

Cp  pressure  coefficient 

F  reduced  vorticity  -  defined  by  Eq.  (5) 

C  perturbation  stream  function  -  defined  by  Eq.  (7) 

H  total  head  -  defined  by  Eq.  (90) 

3G 

SL  spline  derivative  approximation  of  -jj— 

%2q 

L  spline  derivative  approximation  of 

L  body  length  (dimensional) 

number  grid  intervals  in  n-direction 
number  grid  intervals  in  ^-direction 
r  radial  distance 

rB  body  radius 

♦ 

rE  outer  cylinder  radius  (where  far  field  boundary  condition 
applied) 

u  velocity  component  in  x-direction 

v  velocity  component  in  r-direction 

x  axial  distance 

Xq  axial  coordinate  of  initial  line  in  computational  domain 

x^  axial  coordinate  of  final  line  in  computational  domain 

* 

free-stream  speed  (dimensional) 

An  step  size  in  n-direction 

Af  step  size  in  ^-direction 

n  transformed  radial  coordinate 

£  transformed  axial  coordinate 

j;  vorticity  magnitude 


All  other  quantities  are  defined  in  the  text. 
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NOMENCLATURE  (cont.) 

All  quantities  in  the  text  are  made  dimensionless  as  follows: 

* 

distances  by  L 

* 

velocities  by  U 

■k  ^  & 

stream  function  by  L  U 

J  CO 

•k  k 

vorticity  by  U^/L 
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I.  INTRODUCTION 

The  calculation  of  the  flow  field  in  the  vicinity  of  the  tail  of  a 
propeller/body  combination  is  complicated  by  three  effects  which  interact 
strongly  with  one  another.  These  effects  are: 

1.  The  curving  outer  streamlines  adjusting  to  the  upcoming  wake 
line  of  symmetry  induces  a  significant  pressure  gradient  in  the 
boundary  layer  normal  to  the  main  flow  direction. 

2.  The  outer  flow  (where  vorticity  is  negligible)  is  displaced  by 
the  rapidly  thickening  boundary  layer  near  the  tail. 

3.  The  streamlines  near  the  tail  are  modified  by  the  presence  of 
the  propeller. 

The  treatment  of  normal  pressure  gradient  and  displacement  effects 
has  usually  involved  retaining  the  division  of  the  flow  field  into  a 
boundary  layer  and  inviscid  flow.  These  effects  are  then  accounted  for 
approximately  by  the  displacement  body  approach  where  the  boundary-layer 

A 

displacement  thickness  is  added  to  the  original  body  [1  -  3],  or  by 
patching  the  boundary  layer  to  the  inviscid  outer  flow  at  the  edge  of 
the  boundary  layer  [4,  5].  In  either  case  the  solution  must  be  found  by 
iteration  because  of  an  initially  unknown  boundary  location.  The  solution 
process  is  complicated  by  the  iteration  being  notoriously  sensitive  in 
the  tail  region  where  streamline  curvature  is  greatest.  Consequently 
such  a  procedure  requires  a  great  deal  of  judgment  on  the  part  of  the 
user. 

The  effects  of  a  propeller  are  often  accounted  for  by  a  frozen 
vorticity  calculation  in  conjunction  with  an  actuator  disk  to  represent 
the  propeller.  The  boundary  layer  establishes  the  vorticity  distribution 
* 

Numbers  in  brackets  designate  References  at  end  of  text. 
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at  some  station  upstream  of  the  propeller.  Then  the  flow  is  calculated 
through  the  actuator  disk  as  though  it  were  inviscid  but  rotational  - 
hence  the  name  "frozen"  vorticity.  The  result  is  a  modified  flow  field 
which  approximately  accounts  for  the  influence  of  the  propeller.  Huang, 
et  al.  [2]  have  combined  this  type  of  calculation  with  the  displacement 
body  idea. 

The  first  treatment  of  all  three  effects  together  where  simplifying 
assumptions  were  held  to  a  minimum  was  by  Schetz  and  Favin  [6,  7],  Their 
approach  is  based  on  the  full  axisymmetric  time-averaged  unsteady  Navier- 
Stokes  equations  with  the  propeller  modeled  by  an  actuator  disk  or  zone 
and  the  turbulence  field  modeled  by  a  simplified  turbulent  kinetic  energy 
transport  equation  which  depends  only  on  the  streamwise  distance. 

The  solution  strategy  followed  by  Schetz  and  Favin  was  to  solve  the 
Navier-Stokes  equations  only  where  absolutely  necessary,  i.e.  in  a 
restricted  domain,  as  shown  in  Fig.  1,  whose  size  is  dictated  on  physical 
grounds.  Referring  to  Fig.  1,  in  domain  1  the  effects  of  viscosity  are 
negligible  so  that  the  potential  flow  approximation  is  valid,  in  domain  2 
the  conventional  first-order  boundary-layer  equations  hold  and  in  domain  3 
the  full  Navier-Stokes  equations  must  be  used.  Domain  3  is  the  region 
where  the  three  strong  interaction  effects,  mentioned  previously,  are 
present. 

Because  cylindrical  coordinates  are  used,  the  solutions  obtained  by 
Schetz  and  Favin  were  limited  to  bodies  with  tails  which  are  either 
infinitely  thin  or  conical.  The  conical  case  was  accommodated  by  using 
a  stair-step  grid  at  the  body  surface  which  does  not  lead  to  effective 
step  size  control  in  the  boundary-layer  region.  The  equations  of  motion 
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were  solved  in  stream  function-vorticity  form  using  a  direct  solver  for 
the  Poisson  equation  and  the  ADI  technique  for  the  time-dependent  vorticitv 
transport  equation. 

The  aim  of  the  present  work  is  also  to  solve  the  strong  interaction 
problem  associated  with  a  propeller/body  combination  with  as  few  simpli¬ 
fications  as  possible  but  with  more  generality  than  the  approach  of  Schetz 
and  Favin.  The  present  method  is  designed  to  handle  an  arbitrary,  smooth, 
pointed  tail  shape  on  which  is  mounted  either  an  open  propeller  or  a  pump- 
jet  (ducted  propeller) . 

The  solution  strategy  to  be  used  is  very  similar  to  that  of  Murphy  [8] 
in  his  development  of  an  efficient  Navier-Stokes  solver  and  consists  of 
the  following  five  elements: 

1.  A  truncated  computational  domain  (similar  to  Schetz  and  Favin) 
is  used  in  the  tail-near  wake  region  in  which  the  "thin-layer" 
version  of  the  time-averaged  Navier-Stokes  equation  is  solved. 

The  thin-layer  approximation  involves  neglect  of  viscous/turbulent 
diffusion  terms  in  the  streamwise  direction. 

2.  The  steady-state  form  of  the  equations  of  motion  are  solved 
which  usually  involves  many  fewer  iterations  than  obtaining  a 
steady-state  solution  from  the  unsteady  equations. 

3.  The  propeller  is  modeled  by  an  actuator  disk  or  zone. 

4.  A  non-orthogonal  coordinate  transformation  is  used  so  that  an 
arbitrary,  smooth,  pointed  tail  body  and  actuator  disk  or  zone 
are  families  of  the  coordinate  system.  Thus  boundary  conditions 
are  easy  to  apply  and  effective  step  size  control  can  be  main¬ 
tained  close  to  the  body  surface  where  flow  gradients  are  large. 
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5.  A  fourth-order  spline  or  sixth-order  Hermite  discretization  is 
used  in  the  direction  normal  to  the  main  flow  so  that  the  region 
of  high  flow  gradients  near  the  body  surface  can  be  resolved 
with  relatively  few  node  points.  The  result  is  a  more  computa¬ 
tionally  efficient  method.  Murphy  used  a  generalized  Galerkin 
technique  with  splined  Taylor  series  expansions  to  achieve 
fourth-order  accuracy  in  the  direction  normal  to  the  main  flow. 
The  present  approach  achieves  essentially  the  same  result  but 
requires  less  algebra  to  arrive  at  the  numerical  algorithm. 

The  work  is  being  performed  in  two  phases.  In  phase  1,  covered  in 
this  report,  viscosity  and  the  actuator  zone  are  omitted  to  simplify  the 
algorithm.  Otherwise  all  the  elements  of  the  numerical  strategy  are 
present.  The  intent  is  to  gain  an  understanding  of  the  proposed  method 
by  solving  a  simpler  sub-problem  whose  solution  will  then  be  used  as  a 
first  guess  for  the  full  problem.  Thus  in  phase  1,  the  problem  reduces 
to  the  solution  of  Poisson's  equation  for  the  stream  function  with  a 
"frozen"  vorticity  distribution  established  on  the  upstream  computational 
boundary  by  the  boundary-layer  solution.  In  phase  2,  all  of  the  elements 
of  the  problem  will  be  included.  In  this  phase  the  turbulence  field 
will  be  modeled  by  a  simple  two-piece  algebraic  eddy  viscosity  model 


corrected  for  streamline  curvature  effects. 
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II.  ANALYSIS 

The  Frozen  Vorticity  Equation 

In  cylindrical  coordinates  (x,  r,  0)  with  velocity  components  (u,  v,  0) 
where  3/39  E  0  (axial  symmetry),  the  vorticity  vector  has  a  component  in 
the  0-direction  only  with  magnitude  £  given  by 


3v  3u 
3x  3r 


(1) 


For  steady  incompressible  flow  the  equation  of  mass  conservation  reduces 
to 


37  (ru)  +  37  (rv)  =  0  * 


(2) 


Mass  conservation  is  satisfied  identically  by  the  Stokes  stream  function 
defined  by 


3i|;  dtli 

ru  =  37  •  rv  "  -  37 


(3) 


With  the  aid  of  the  Stokes  stream  function  the  vorticity,  Eq.  (1) , 
becomes 


_3jJ)  1  dijj  3 _ 

3r2  "  r  9r  3x2 


(4) 


For  axially  symmetric  inviscid  rotational  flow  it  can  be  shown  that  [9] 


r,  =  r  F(< P) 


(5) 


m 
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where  F  (ip)  is  called  the  reduced  vorticlty  and  is  determined  by  conditions 
at  the  upstream  boundary.  In  terms  of  F(ip),  Eq.  (4)  becomes 


1  9ijj  djp 

r  3r  „  2 

dx 


r2  F(ip) 


(6) 


In  the  far  field  the  flow  must  approach  a  uniform  stream,  i.e.  u~>-l 

2 

as  r-*50.  Thus  the  stream  function  must  behave  like  r  /2  as  r  goes  to 
infinity.  If  this  unbounded  behavior  is  subtracted  from  ip  the  result  is 
a  perturbation  stream  function  G(x,  r)  defined  by 

2 

G(x,  r)  =  ip(x,  r)  -  —  ,  (7) 

2 

which  is  well  behaved  as  r-*».  Since  r  /2  is  a  particular  solution  of 
Eq.  (4),  G  satisfies  the  same  differential  equation  as  ip,  viz.: 


1  9G  d2G 

r  9r  _  2 

dx 


r2  FOP) 


(8) 


For  the  truncated  computational  domain  shown  in  Fig.  2  the  appropriate 
boundary  conditions  are: 


4)  =  0  on  r 


10 


X0  -  x  -  1 


x  >  1 


(9) 


ip  =  \pn(r)  on  x  -  xQ  ,  rB  1  r  1  rE 


(10) 


x^xo 


^  on  r  =  rE 


(ID 
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The  boundary  condition  on  x  =  x^,  the  right-hand  boundary,  is  somewhat 
arbitrary  and  could  just  as  well  be  replaced  by  \pxx  =  0. 

Body  Fitted  Coordinates 

For  ease  of  application  of  boundary  conditions  a  coordinate  transfor¬ 
mation  must  be  found  which  maps  the  body  surface  and  propeller  shroud 
into  coordinate  surfaces  in  the  mapped  plane.  In  addition,  for  modeling 
the  propeller  by  an  actuator  disk  or  zone,  vertical  lines  should  remain 
unchanged.  The  following  sheared  transformation  has  these  desired 
properties : 


£  =  x  (13a) 

n  =  V50  '  r  -^2(x)  (13b) 

where  9^(x)  an<i  (?2^x^  are  piecewise  continuous  functions  which  depend  upon 

the  particular  geometry.  The  transformation  applies  to  a  shroud  with  or 

without  thickness  whose  leading  and  trailing  edges  are  sharp,  i.e.,  which 

has  an  included  angle  considerably  less  than  90  degrees.  For  blunt  leading 

or  trailing  edges  the  transformation  breaks  down. 

Figure  3  illustrates  a  typical  application  of  F.q.  (13)  to  a  body- 

shroud  combination.  Because  of  the  particular  geometry  the  region  must 

be  divided  into  five  parts.  As  can  be  seen,  the  body  is  mapped  onto  part 

of  n  =  0  and  the  shroud  becomes  a  slit  on  n  =  n  .  Figure  3  also  illus- 

s 

trates  that  lines  of  constant  £  and  constant  n  in  the  physical  plane  are 
generally  not  orthogonal. 

Since  the  derivatives  of  <().^(x)  and  ^(x)  are  discontinous  at  region 
boundaries,  the  transformed  differential  equation  will  have  coefficients 
which  are  also  discontinuous  at  region  boundaries.  Continuity  of 

L  . 
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the  solution  iJj(£.  0)  and  its  derivatives  will  provide  the  necessary 
"jump"  conditions  at  the  boundary,  as  will  be  seen. 

In  the  present  "test"  problem  we  restrict  the  geometry  to  a  body  with 
no  shroud  so  that  only  two  mapped  regions  are  necessary,  as  illustrated 
in  Fig.  2.  The  appropriate  transformations  are: 

Region  I:  x^  <_  x  1  (aft  portion  of  body) 

£  =  x  ,  (14a) 


n 


r 


rR(x) 

-V 


(14b) 


Region  II:  x  >  1 


£  =  x 


(15a) 


n 


(15b) 


In  region  II  we  are  merely  scaling  r  to  be  compatible  with  region  I. 

The  transformation  in  II  is  therefore  orthogonal.  In  region  I  we  see 
that  n  =  0  on  r  =  r  while  in  region  II  ri  =  0  on  r  =  0.  In  both  regions 

D 

r|  =  1  on  r  =  r_. 

E 

Transformed  Equation  of  Motion 

By  applying  the  chain  rule,  F.q.  (8)  in  transformed  (£,  r|)  coordinates 
can  be  written  in  the  following  general  form  valid  for  regions  I  and  II: 
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Region  II: 

A  =  Ar  =  =  0  ,  (28) 

B  =7-  ,  (29) 

F. 

Bn  =  0  .  (30) 

From  Eqs.  (9)  -  (12)  together  with  Eq.  (7),  the  boundary  conditions 
for  the  transformed  equation  of  motion  are: 


G  =  “  \  rj(.r,)  ,  on  n  =  0  ,  c,Q  <  T,  <_  1  ,  (31) 

C  =  0  ,  on  ri  =  0  ,  (  >  1  ,  (32) 

0  =  8(H)  ,  on  ?  =  ro  ,  0  <_  n  <  1  ,  (33) 

c  =  CE(0  ,  on  n  =  1  ,  X  r0  ,  (34) 

~  =*  0  ,  on  K  =  F  .  0  1  n  1  1  ■  (35) 


The  boundary  conditions  must  be  supplemented  by  "jump"  conditions 
on  the  derivatives  of  G  in  the  transformed  plane  at  F  =  1  so  that  the 
solution  can  be  continued  from  one  region  to  the  next.  We  know  that  G 
and  its  x  and  r  derivatives  must  be  continuous  at  C  =  1  because  the 
stream  function  and  velocity  components  are  continuous.  Thus  at  £  =  1 


we  have 
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(G)  =  (G) 


(G  )  =  (C  1 

x  —  x  + 


«V-  ‘  <Cr>+ 


at  x  =  1 


(36) 

(37) 

(38) 


where  denotes  the  left-hand  limit  (in  region  I),  "+"  denotes  the 
right-hand  limit  (in  region  II)  and  the  tilde  denotes  the  function  in 
region  II.  The  chain  rule  relating  first  derivatives  in  (£,  n)  variables 
to  those  in  (x,  r)  variables  can  be  expressed  as 


+  A(^,n) 


(39) 


B(<r,n) 


(40) 


Then  with  the  aid  of  the  chain  rule  relations,  Eqs.  (39)  -  (40),  and  the 

mapping  derivatives,  Eqs.  (23)  -  (24)  and  (28)  -  (29),  and  noting  that 

B  (l,n)  =  B  (l,n)  since  rD(l)  =  0,  the  jump  conditions  can  be  written 
X  IX  D 

as 


(Gf)+  -  (C  )_  +  Ajd,  n)(on>_  ,  (41) 

■  (V-  •  <42) 


From  Eq . 


(42)  and  the  fact  that  B^ 


=  0  in  regions  I  and  II  we  also  have 
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We  can  infer  in  this  case  that  all  the  n-der ivatives  are  continuous  at 
£  =  1.  Thus  the  only  "jump"  occurs  in  the  ^-derivative,  as  given  by 
Eq.  (41). 

Numerical  Algorithm 

For  the  frozen  vorticity  equation  the  numerical  algorithm  has  been 
chosen  to  be  almost  the  same  as  the  one  to  be  used  later  on  for  the  full 
problem  with  viscous/turbulent  effects.  The  reason  for  this  choice  was 
to  gain  familiarity  with  the  algorithm  on  a  simpler  problem. 

A  combination  finite  difference  and  spline  representation  is  used 
for  reasons  of  accuracy  and  maintenance  of  a  block  tridiagonal  matrix 
structure  of  the  resulting  equations.  In  the  streamwise  direction, 
where  flow  quantities  are  expected  to  vary  somewhat  slowly,  centered 
difference  formulas  with  constant  A£  are  used.  In  the  direction  normal 
to  the  £-axis,  where  flow  quantities  are  expected  to  vary  rapidly  near 
the  body  (in  the  viscous/turbulent  case),  polynomial  spline  formulas 
are  used  to  obtain  higher-order  accuracy. 

We  now  define  the  following  spline  derivatives  in  the  n-direction: 

(44) 

Then  the  frozen  vorticity  equation  can  be  written  at  node  point  (i,  j)  as 
<“  L  +  b  \  +  Ca  ♦  d  l  +  e)Jf  }  -  0  ,  <«> 


where  i  and  j  denote  the  £  and  n  coordinate  lines  respectively.  The 
centered  difference  approximations  in  the  f-direction  are: 
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(47) 


G,  ,  .-20  .  +  G...  . 

(r  )  c  — L j _ _ L’ J _ 1±L — 1  + 

1»  1  a r 


o(A  r2) 


(48) 


Upon  substitution  of  the  central  difference  (CD)  approximations  in  Eq .  (46) 
an  algebraic  expression  results  for  0,  £,  and  L  at  the  nodal  points 
(i-1,  j),  (i,  j),  and  i+1,  j).  For  the  method  of  solution  chosen,  spline 
line  overrelaxation  (SPI.0R),  experience  has  shown  that  sweeping  in  the 
direction  is  much  more  straightforward  than  in  the  n -direction  because 
of  the  piecewise  continuous  mapping. 

For  a  "-sweep  using  SPLOR  the  unknowns  lie  along  line  i  with  those 
on  line  (i-1)  and  (i+1)  considered  known.  Thus  from  Eq .  (46),  together 
with  the  CD  approximations,  the  unknowns  at  line  i  are  found  to  be: 


-  G.  .  +  d4  .  £.  .  +  a,  .  L  .  =  P.  . 

2  l>j  i.  .i  ij.i  i»j  i».i  i>.i 

where. 


(49) 


bi  1 

P  =  — =-* — L  _  o  )  - 

ri,  j  2 A ■’  1  i-1,  j  i+1 , 


Ac 


2  j  +C1+1.  J> 


a,  ,  .  (50) 
1  »  .1 


Because  splines  are  being  used  in  the  n-direction  Eq .  (49)  involves 
unknowns  only  at  j.  Equation  (49)  is  valid  in  both  regions  1  and  II 
except  at  n  =  0  where  the  coefficient  d  becomes  unbounded.  At  n  =  0  a 
limiting  form  of  Eq.  (49)  must  be  found.  In  addition,  a  special  form 
of  Eq.  (49)  must  be  derived  for  the  junction  line  between  regions  I  and 
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The  Junction  Line 

The  junction  line  is  treated  by  a  generalized  version  of  the  method 
originally  developed  by  Chmielewsk.i  and  Hoffman  [10]  for  SLOR  solution  of 
elliptic  partial  differential  equations  with  discontinuous  coefficients. 

In  this  method  the  solution  in  each  region  is  extended  by  analytic  continu' 
ation  one  step  beyond  the  boundary  into  the  other  region.  As  a  result, 
two  lines  of  "fictitious"  unknowns  are  formed,  one  on  either  side  of  the 
junction.  In  the  present  case  the  fictitious  unknowns  can  be  eliminated 
beforehand  by  use  of  the  jump  conditions  and  by  approximating  in 
region  I  by  a  second-order  backward  difference  formula.  As  a  result,  a 
single  algebraic  equation  is  obtained  on  the  junction  line  so  that  SPLOR 
can  still  be  used.  In  the  Chmielewski-Hoffman  method  the  system  of 
unknotfns  on  the  junction  line  and  two  adjacent  fictitious  lines  was  not 
reduced  to  a  single  tri-diagonal  equation.  Consequently  point  relaxation 
was  used  to  solve  this  system,  but  everywhere  else  SLOR  was  used.  We 
note  that  by  proper  matrix  partitioning  the  C-H  method  becomes  block 
tridiagonal  so  that  SPLOR  can  be  used  at  the  junction  line. 

On  the  junction  line  between  regions  I  and  II,  denoted  by  i  =  IJ, 


a 


5 


in  region  I  is  approximated  by  a  second-order  backward  difference 


formula,  viz.: 


+  o(Af;2) 


(51) 


A  centered  difference  formula  at  i  =  IJ  is  still  used  for  G^r.  Thus  in 
region  I  the  transformed  equation  of  motion,  Eq.  (46),  becomes  at  (IJ,  j): 
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a  ,  ,  ,  L*  •  .  bx  t  • 

IJ,  j  IJ,  j  IJ,  J 


iT.  0  .  -  4  JtTT  .  .  +  3  £tt  7 

IJ-2,  1 _ IJ-1.  3 _ IJ,  1 


2A£ 


CIJ-l  -  2  CIJ ,  J  +  W.  j  +  d  ,  .  +.  .-0 

Ar2  IJ,  J  IJ,  .1  IJ,  J 


where  the  asterisk  denotes  a  fictitious  unknown  from  the  analytic  continu¬ 
ation.  Solving  for  unknowns  on  i  =  IJ  and  IJ  +  1,  we  get 


£>-2  GIJ,  j  +  ^2  GIJ+1,  j  + 


u,  j  2AC 


IJ,  j 


+  aIJ,  j  LIJ,  j  °j 


(52) 


where  Q.  consists  of  quantities  considered  known  and  is  given  by 


Qj  " 


IJ,  j 


\r 2  U-l,  j 


IJ,  i 
2A  F, 


(£U-2,  j  -  4 


eU-l, 


(53) 


Similarly,  in  region  II  Eq.  (46)  yields  the  following  relation  for 
unknowns  on  i  =  IJ  and  the  fictitious  unknowns  on  the  line  i  =  IJ  -  1: 


A? 

where 


1  j.  ~ 

2  GIJ,  j  +  ^2  GIJ-1,  j  +  dIJ,  j  £IJ,  j  +  3IJ,  j  LIJ,  j 


(54) 


Qj  "  GIJ,  j 


Af 


2  IJ+1,  j 


(55) 


In  obtaining  Eq.  (54)  we  have  made  use  of  the  fact  that  b  =  0  in  region  II. 

To  obtain  a  single  equation  at  (IJ,  j)  we  make  use  of  the  continuity 
and  jump  conditions  at  i  =  IJ  to  eliminate  unknowns  which  arise  from  the 
analytic  continuation.  The  continuity  conditions  were  found  to  be 
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GIJ,  j  =  f’U,  j 
?IJ,  j  =  £IJ,  j 
LIJ,  j  =  LIJ,  j 


and  the  jump  condition  was 


(Gr)TT  .  =  (Gr)TT  .  +  Att  . 
C  Id  »  j  S  Id  ,  j  Id  ,  J 


“Id,  j 


Upon  replacement  of  the  C-derivatives  by  centered  differences,  the  jump 
condition  becomes 


Gid+i,  i  ~  Gid-1,  j 

2AC 


Gld+1.  i  ~  Gld-1,  i 

2AC 


j  Id,  j 


(56) 


By  combining  the  continuity  conditions  and  jump  condition,  as  given  by 

~  -k 

Eq.  (56),  with  Eqs .  (52)  and  (54)  we  can  eliminate  the  unknowns  GIJ_1  ^ 
and  G_  _  .  .  and  arrive  at  the  following  equation  involving  only 

J 

unknowns  on  i  =  I J : 


AC 


-  Gtt  .  +  d.  S.tt  .  +  a.  Lt  .  =  P. 
2  IJ,  j  J  IJ,  J  J  Id,  J  .1 


where: 


Sj  '  2  <aIJ,  j  *  5IJ,  j>  ' 


3,  '  1  <du,  j  +  d»,  +  -fe1 


Pi  '  eIJ,  j  ‘  Af2  j  +  CU+] ,  j’ 


IJ,  d 


a. 


4AC  v~IJ-2,  j  4  *td-l,  j} 


(57) 

(58) 

(59) 


(60) 
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Equation  (57)  is  seen  to  be  of  the  same  form  as  Eq .  (49)  for  the  general 
case  with  a  modification  of  the  coefficients  due  to  the  discontinuity  in 
the  mapping  derivatives. 

Spline  Formulas 

With  the  introduction  of  the  spline  derivatives  Z  and  L  we  have  seen 
that  the  frozen  vorticity  equation  yields  a  single  equation  for  the 
unknowns  G,  l,  and  L  at  a  node  point  (i,  j).  To  complete  the  set  of 
equations  at  (i,  j )  two  polynomial  spline  formulas  relating  the  function 
and  its  n-derivatives  must  be  introduced.  For  the  sake  of  simplicity 
we  will  use  spline  relations  of  fourth-order  accuracy  with  variable  Aq 
as  given  by  Rubin  and  Khosla  [11].  As  a  further  simplification  we  will 
use  the  same  spline  formula  twice,  namely  5^(4,  0)  (in  the  nomenclature 
of  Rubin  and  Khosla)  to  relate  Z  to  G  and  L  to  Z,  viz.: 


and 


3.  G.  .  .  +  a.  G.  .  +  y.  G.  ... 
J  i,  j-1  J  i,  J  J  i,  1+1 


+  a2  ZJ  .  ,  +  (1  +  a)2  Z.  .  +  Z.  =  0 

i»  1-1  i»  J  i»  J+l 


3.  Z.  .  +  a.  Z  .  +  y.  Z. 

J  i,  j-1  J  i,  J  J  ii  J+l 


+  a2  L,  .  .  +  (1  +  a)2  L,  .  +  L.  =  0 

i.  J-1  1 >  J  i.  J+l 


(61) 


(62) 


where  a,  3,  Y>  and  a  are  functions  only  of  the  step  size  and  are  given 

by 
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a. 

.1 


2(1  -  a)  (1  +  n',? 

An.  ,  o 

.1-1 


p  =  2ct?(2  +  o) 

1  An.  ,  O  +  a) 

1-1 

_  2(1  +  2a) 

Yi  An .  , (1  +  0)0 
1-1 


(63) 


(64) 


(65) 


and 


0 


0. 

1 


AVl 


(66) 


An. 

i 


'.1+1 


n. 

i 


(67) 


Other  polynomial  soline  approximations  can  be  used  in  place  of  F.qs.  (61) 
and  (62)  depending  upon  the  accuracy  desired. 

In  the  (i,  j)  grid  point  counting  scheme  adopted  here,  i  =  1  denotes 
£  =  and  i  =  M  +  1  denotes  t,  =  while  1  =  1  denotes  0  =  0  and  j  =  N  +  1 
denotes  n  =  1.  Thus  the  number  of  grid  cells  is  M  by  N. 


Conditions  on  Lower  and  Upper  Boundaries 

At  each  boundary  (n  =  0  and  1)  three  equations  are  needed  to  determine 
the  unknowns  G,  i,  and  L.  The  first  equation  is  derived  from  the  boundary 
condition,  the  second  equation  from  the  equation  of  motion  and  the  third 
equation  from  an  appropriate  spline  relation. 

On  n  =  0  the  boundary  condition  can  be  expressed  as 


f(V 


(68) 


where, 


f  I  2  , 

2  rB  ( 


T,  <  1 


=  \ 


(69) 
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In  region  I  the  second  relation  on  0  =  0  is  provided  by  Eq .  (49) 
evaluated  at  j  =  1,  viz. 


— — r  G.  .  +  d.  ,  Z.  1  +  a.  L .  =  P.  1  ,  <  ?  <  1 

.-2  i,l  i,  1  l,  1  l,  1  l,  1  l,  1  0  — 


(70) 


In  region  II,  as  well  as  on  the  junction  line,  Eq.  (70)  is  not  valid 
because  d.  -.is  unbounded,  as  already  mentioned.  We  must  therefore 
replace  Eq.  (70)  by  an  appropriate  limit  derived  from  the  equation  of 
motion.  On  q  =  0  the  governing  differential  equation  reduces  to 


L  -  lim 

n 


im  Hr 

— o  W| 


which  can  only  be  satisfied  if  1  =  0.  Thus  the  condition 


*i,  1  =  °  ’  ’ 


(71) 


replaces  Eq .  (70)  in  region  II.  The  final  equation  is  supplied  by  the 
spline  interpolation  polynomial  S^(4,  0),  as  discussed  by  Rubin  and 
Khosla  [11].  The  formula  is 


An 


_  n  _| _ [L-  Q  _  -b—  1  _  .... ■?_  £ 

2  An2  °i,  2  Anx  i,  l  An1  i. 


2  -  L.  1  =  0(An1)  .  (72) 


The  use  of  Eq .  (72)  means  that  G  is  determined  to  fourth-order  accuracy, 
H  to  third-order  accuracy  and  L  to  second-order  accuracy. 

At  the  outer  boundary,  n  =  1,  we  have  the  boundary  condition  common 
to  both  regions  I  and  TI: 


Gi,  N+l  CE. 

i 


(73) 
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The  second  equation  is  the  equation  of  motion  evaluated  at  n  =  1,  viz.: 


A£ 


2  Gi,  N+l  +  di,  N+l  ^i,  N+l  +  ai,  N+l  Li ,  N+l  Pi,  N+l  ’  (74) 


and  the  third  equation  is  a  spline  relation  derived  from  S  (4,  0) , 
similar  to  Eq .  (72): 


a. 


+  L. 


2  i,  N  An  i,  N  2  ui,  N+l  An  *1,  N+l  “i,  N+l 
N  mN 


0  .  (75) 


Matrix  Equations 

A  compact  form  of  the  spline-finite  difference  equations  is  obtained 
by  introducing  the  following  three-component  column  vector  of  unknowns: 


r  o  i 


z.  . 

1»  J 


9. 


L  I.  . 
-  i»  J 


(76) 


In  terms  of  Z.  .  the  equations  in  the  field  and  on  the  boundaries 
1  j  J 

associated  with  line  i  can  he  written  as  a  block  tridiagonal  system,  viz. 


A1  Z1  +  C1  Z2 


B.  Z.  .  +  A.  Z.  +  C.  Z  .  .  =  R . 
J  .1-1  J  J  J  J+l  .1 


bn+i  zn  +  an+i  zn+i  rn+i 


2  <  j  <  N 


(77) 


where  the  i  index  is  understood.  In  Eq.  (77)  the  coefficients  A,  B,  C 
are  3x3  matrices  which  depend  only  on  body  geometry  and  step  size,  and 
Rj  is  a  three-component  column  vector  of  known  quantities  which  is 
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obtained  from  the  various  field  and  boundary  equations.  The  linear 
system,  Eqs.  (77),  can  be  solved  by  standard  lower-upper  (L-U)  factorization 
methods  -  see  Keller  [12]. 

Initial  and  Final  Lines 

To  start  the  calculation,  data  on  the  initial  line,  i  =  1,  are 
required,  namely:  G,  £,  and  L  at  each  node  point.  In  addition,  the 
reduced  vorticity  £  and  total  head  H  are  needed.  With  ip  and  u  as  input 
on  the  initial  line,  then  G  and  £  are  computed  from: 


2 

G  =  ’p  -  — j  ,  (78) 

and 

£  =  r(rE  -  rB)(u  -  1)  ,  (79) 

where  Eq.  (79)  is  derived  from  the  Stokes  stream  function  relation 
<p  =  ru.  The  second  spline  derivative  L  on  the  initial  line  is  computed 
by  applying  (4 ,  0)  to  £.  Finally,  the  reduced  vorticity  F  is  determined 
from  Eq.  (5)  and  is  stored  in  a  table  versus  ip  as  is  the  total  head  H. 

At  the  final  line  of  the  computational  domain,  i  =  M  +  1,  the 
boundary  condition  is  tp^  =  0.  Using  a  CD  approximation  this  condition 
becomes 


V?,  j  V  j  ’ 


(80) 


and  hence  becomes  (note  that  in  region  II,  b  =  0) : 


M+l,  j 


.1  UM,  j  "  CM+1,  j 


A.c 


(81) 


Other  than  this  modification  the  calculation  at  i  =  M  +  1  is  the  same  as 


at  any  other  line. 
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III.  RESULTS 

Potential  Flow  Solutions 

As  a  means  of  testing  the  present  method,  potential  flow  solutions 
were  generated  and  compared  with  results  obtained  from  the  second-order 
surface  singularity  method  of  Fernandez  [13].  Two  body  shapes  were  run: 
the  modified  spheroid  of  Patel  [14]  and  the  F-57  low-drag  body  of  Parsons 
and  Goodson  [15]. 

The  SPLOR  solutions  for  the  two  test  cases  were  obtained  with  grid 
parameters  as  given  in  Table  1. 


Modified  Spheroid 

F-57  Body 

xo 

0.4816 

0.2280 

X1 

2.0368 

1.7720 

rE 

3.4215 

3.2430 

N 

n 

20 

20 

Nc 

60 

40 

Anl 

0.0054 

0.0057 

cn 

1.2 

1.2 

Table  1.  Potential  Flow  SPLOR  Grid  Parameters. 

To  simplify  the  input  G^(^)  was  set  to  a  constant  value  and  r^  taken 
large  enough  so  that  this  approximation  had  negligible  effect  on  the 
solution.  A  highly  nonuniform  step  size  in  u  was  used  to  test  this 
feature  of  the  algorithm.  As  shown  in  Table  1,  the  initial  An 
(adjacent  to  the  body)  was  very  small  with  the  ratio  of  successive  steps 


in  n  held  fixed,  i.e., 
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An. 


4nj-i 


constant 


An  overrelaxation  factor,  r^,  of  1.6  was  used  throughout  the  calculations 

on  G,  £,  and  L  with  no  attempt  made  to  determine  the  optimum  r^.  Runs 

were  also  made  with  smaller  step  sizes  than  those  shown  in  Table  1  but 

the  change  on  the  numerical  solution  was  negligible. 

The  numerical  solution  was  considered  converged  when  the  largest 

residual  change  in  G,  i.e.  |<5G^|  =  |G^N+^  -  G^  |  ,  was  less  than  a 

prescribed  amount,  usually  5  x  10  In  determining  convergence  the 

largest  residual  change  in  SL  or  L  could  just  as  well  have  been 

monitored.  As  would  be  expected,  for  the  same  number  of  steps  in  r, 

those  runs  with  a  highly  nonuniform  distribution  took  more  iterations 

than  those  with  a  uniform  distribution  to  arrive  at  the  same  cutoff 

tolerance.  Typical  number  of  iterations  to  reach  a  tolerance  on  'G 
_  ^ 

of  5  x  10  was  from  40-50  for  the  parameters  given  in  Table  1  with 
CPU  time  on  an  IBM  370/3033  of  about  90  seconds  (using  double  precision 
arithmetic) . 

Figures  4  and  5  show  the  body  pressure  distribution  (x  <  1)  as 

well  as  the  centerline  pressure  distribution  downstream  of  the  body 

(x  >  1)  for  the  modified  spheroid  and  F-57  respectively.  The  surface 

singularity  solution  was  obtained  only  on  the  body.  Comparison  of  SPLOR 

with  surface  singularity  results  shows  excellent  agreement  except  very 

close  to  the  tail.  The  SPLOR  value  of  C  at  the  tail  point  is  seen  to 

P 

be  too  low.  Better  agreement  with  the  surface  singularity  results  by 
SPLOR  would  probably  be  obtained  by  a  nonuniform  grid  in  ^  which  is  denser 
near  the  tail  point  but  this  is  not  allowed  in  the  present  formulation. 


L 
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Frozen  Vorticity  Test  Case 

The  F-57  low-drag  body  was  used  as  a  test  case  for  a  frozen  vorticity 
calculation.  The  initial  profile  was  chosen  at  x^  =  0.7684  which  is  at 
the  approximate  beginning  of  the  strong  interaction  region. 

The  thickness  of  the  vortical  layer  was  taken  to  be  -  0.12  which 
corresponds  approximately  to  the  thickness  of  the  turbulent  boundary 
layer  at  the  initial  station  for  a  Reynolds  number  based  on  body  length 
of  106. 

The  computed  turbulent  boundary-layer  profile  at  x^  =  0.7684  was  not 
used  for  the  initial  frozen  vorticity  profile  because  the  values  of 
vorticity  near  the  wall  were  considered  to  be  too  large  and  the  no-slip 
condition  at  the  wall  is  inappropriate.  Instead,  a  "manufactured" 
velocity  profile  was  generated  which  resembled  the  turbulent  boundary- 
layer  profile  but  had  lower  values  of  vorticity  near  the  wall  and  a 
wall  slip  velocity. 

Some  care  had  to  be  exercised  in  order  to  make  the  manufactured 
velocity  profile  compatible  with  the  potential  flow  at  the  edge  of  the 
vortical  layer  (n  =  n^) .  The  procedure  finally  found  to  be  satisfactory 
is  as  follows: 


1.  We  assume  that  £  varies  parabolically  with  r| : 

2 

£  =  a±  +  a2n  +  a3n 

where  u  is  related  to  £  by 


u 


1  + 


r(rE  - 


V 


(82) 


(83) 
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and  the  coefficients  a^,  a^,  and  a^  are  determined  from  the  following 
conditions: 


1(0 ) 
i(n5) 

*'(n6) 


'c.q  ,  given  wall  slip  velocity. 

2^  ,  u  must  match  with  the  potential  flow  value  at 

n  =  ns. 

,  the  slopes  of  the  vortical  and  potential  flow 
profiles  must  be  the  same  at  n  *  n^. 


The  resulting  initial  profile  for  u  is  shown  in  Fig.  6. 

2.  In  the  vortical  layer  G  is  then  obtained  by  integration  of  £,  which 
gives 


G  =  ~  \  rB2  +  ain  +  \  a2T]2  +  I  a3q3  *  (84) 

3.  The  potential  flow  must  be  adjusted  for  f|  >  to  account  for  the 
vortical  layer  displacement  effect.  This  correction,  which  must  be 
applied  at  each  nodal  point  in  the  potential  flow  portion  of  the  profile, 
is 


G  ,  =  G  '  .  -  +  AG  , 

corrected  potential  o 


(85) 


where 


AG,  =  (G  _  -  G  _  ,) 

o  vortical  potential  n  =  n. 


(86) 


4.  We  then  obtain  the  vorticity  distribution  by  making  the  boundary- 


layer  assumption. 
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Making  use  of  Eqs.  (14b)  and  (83)  the  vorticity  £  is  found  to  be  related 
to  i  and  L  by 


1 

f  L 

r(rE  "  rB)  1 

rE  ~  rB 

(87) 


where  L  is  obtained  by  differentiation  of  Eq .  (82),  viz.: 

L  =  |^  =  a2  +  2a3n  . 


(88) 


The  vorticity  £  is  generally  not  zero  at  the  edge  of  the  vortical  layer 
because  of  the  neglect  of  3v/9x  in  the  definition.  To  adjust  the  vorticity 
to  zero  at  n  =  a  linear  correction  with  n  is  applied.  The  corrected 
vorticity  C  is  therefore: 


e  -  c  -  ?<n*)  *  — 


(89) 


5.  The  total  head  distribution  in  the  vortical  layer  is  determined  from 
the  dimensionless  Bernoulli  equation: 


2  2 

H  =  C  +  u  +  v 
P 


(90) 


In  the  vortical  layer  H  <  1  whereas  in  the  potential  flow  H  =  1.  In  the 
vortical  layer  u  is  known  from  Eq ,  (83)  but  and  v  are  known  only  at 
the  wall  and  at  n  =  •  To  obtain  and  v  within  the  vortical  layer  we 

assume  that  each  varies  linearly  with  n  which  then  allows  the  initial 
values  of  H  to  be  computed. 
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A  highly  nonuniform  grid  spacing  in  n  was  used  with  =  26  and 
-4 

An^  »  1.5  x  10  .  This  distribution  gave  nine  points  in  the  vortical 

layer. 

A  comparison  of  the  potential  and  frozen  vorticity  pressure 
distribution  on  the  body  and  wake  centerline  as  obtained  by  the  SPI.OR 
method  is  shown  in  Fig.  7.  Also  shown  are  the  experimental  data  of 
Patel  and  Lee  [16]  at  a  Reynolds  number  of  1.2  x  10^.  The  frozen  vorticity 
solution  generally  follows  the  experimental  trend  as  the  tail  is  approached 
and  has  a  maximum  Lower  than  that  of  the  potential  flow  solution. 
Downstream  of  the  tail  the  frozen  vorticity  C^  decays  too  slowly  due  to 
the  presence  of  a  frozen  vorticity  "pressure  wake"  illustrated  by  the 
pressure  profile  plot  at  x^  =  1.309  in  Fig.  8.  The  pressure  wake  arises 
as  a  result  of  the  total  head  defect  on  the  initial  calculation  surface 
together  with  the  requirement  that  H  remain  constant  on  streamsur faces . 

We  empnasize  that  the  point  of  the  present  test  case  is  to  demon¬ 
strate  that  this  portion  of  the  algorithm  is  working  properly.  Conse¬ 
quently  no  attempt  has  been  made  to  explore  the  characteristics  of 
boat  tail  frozen  vorticity  solutions. 
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IV.  CONCLUSIONS 

The  present  application  of  spline  line  overrelaxation  (SPLOR)  in 
conjunction  with  a  piecewise  continuous  body  fitted  coordinate  system  is 
the  first  step  is  applying  this  method  to  the  calculation  of  incompressible 
axisymmetric  viscous/ turbulent  flow  fields.  The  generally  good  agreement 
of  present  results  with  those  of  a  second-order  surface  singularity  method 
for  the  two  potential  flow  cases  presented  gives  some  confidence  in  SPLOR. 

An  important  feature  of  the  algorithm  is  that  line  overrelaxation 
is  maintained  at  junctions  in  the  map  where  the  mapping  derivatives 
are  discontinuous.  This  feature  is  a  result  of  using  backward  differencing 
for  ir.  Another  important  and  well-known  feature  of  the  spline  formulation 
is  that  the  coefficient  matrix  associated  with  each  line  of  unknowns  is 
tridiagonal  which  leads  to  an  efficient  solution  of  the  matrix  equation. 

Experience  with  the  algorithm  in  its  present  form  points  to  several 
shortcomings  which  can  easily  be  corrected.  These  are: 

1.  Use  of  a  constant  step  size  in  the  ^-direction  does  not  allow  for 
proper  resolution  in  the  vicinity  of  the  rear  stagnation  point  unless 
a  very  large  number  of  Intervals  in  £  is  used.  This  deficiency  can 
be  corrected  by  introduction  of  a  stretching  function  in  the  t- 
direction  to  produce  a  nonuniform  grid  near  the  tail,  or  wherever 
else  it  is  needed,  so  that  the  required  accuracy  can  be  achieved 
with  a  minimum  number  of  intervals. 

2.  Application  of  the  far-field  boundary  condition  on  a  cylinder  of 
finite  radius  is  both  cumbersome  and  restrictive.  This  feature  can 
be  removed  by  introducing  an  additional  transformation  in  the  n- 


direction  to  map  an  infinite  radius  into  unity. 


-34-  14  May  1980 

CHH/pjk 


REFERENCES 


1.  D.  F.  Myring,  "The  Profile  Drag  of  Bodies  of  Revolution  in  Subsonic 
Axisvmmetric  Flow,"  Roval  Aircraft  Establishment  Technical  Report 
72234  (1972). 

2.  T.  T.  Huang,  H.  T.  Want,  N.  Santelli,  and  N.  C.  Groves,  "Propeller/ 
Stern/Boundary-Laver  Interaction  on  Axisvmmetric  Bodies:  Theory  and 
Experiment,"  David  W.  Taylor  Naval  Ship  Research  and  Development 
Center  Report  76-0113,  December  1976. 

3.  G.  H.  Hoffman,  "A  Method  for  Calculating  the  Flow  Field  in  the  Tail 
Region  of  a  Body  of  Revolution,"  Journal  of  Ship  Research  (to  appear). 

4.  A  Nakayama,  V.  C.  Patel,  and  L.  I.andweber,  "Flow  Interaction  Near 
the  Tail  of  a  Body  of  Revolution.  Part  2:  Iterative  Solution  for 
Flow  Within  and  Exterior  to  Boundary  Layer  and  Wake,"  Trans.  ASME, 

Series  I,  J.  Fluids  Engin.  £8,  538  (1976). 

5.  H.  E.  H.  Mahgoub  and  P.  Bradshaw,  "Calculation  of  Turbulent  -  Inviscid 
Flow  Interactions  with  Large  Normal  Pressure  Gradients,"  AIAA  Journal, 

17,  1025  (1979). 

6.  J.  A.  Schetz  and  S.  Favin,  "Numerical  Solution  for  the  Near  Wake  of 
a  Body  with  Propeller,"  Journal  Hydronautics  11 ,  136  (1977). 

7.  J.  A.  Schetz  and  S.  Favin,  "Numerical  Solution  of  a  Body  -  Propeller 
Combination  Flow  Including  Swirl  and  Comparisons  with  Data,"  Journal 
Hydronautics  1_3,  46  (1979). 

8.  J.  D.  Murphy,  "An  Efficient  Numerical  Method  for  the  Solution  of  the 
Incompressible  Navier-Stokes  Equations,"  AIAA  paper  77-171,  AIAA  15th 
Aerospace  Sciences  Meeting,  Los  Angeles,  CA,  January  24-26,  1977. 

9.  G.  K.  Batchelor,  "On  Steady  Laminar  Flow  with  Closed  Streamlines  at 
Large  Reynolds  Numbers,"  J.  Fluid  Mech.  ^L,  177  (1  956). 

10.  G.  E.  Chmielewski  and  G.  H.  Hoffman,  "Finite-Difference  Solution  of 

an  Elliptic  Partial  Differential  Equation  with  Discontinuous  Coefficients, 
International  J.  Numer.  Meth.  Engin.  L2,  1407  (1978). 

11.  S.  G.  Rubin  and  P.  K.  Khosla,  "Polynomial  Interpolation  for  Viscous 
Flow  Calculations,"  J.  Comp.  Phvs.  j24,  217  (1977). 

12.  H.  B.  Keller,  "Accurate  Difference  Methods  for  Nonlinear  Two-Point 
Boundary  Value  Problems,"  STAM  J.  Numer.  Anal.  H,  305  (1974). 

13.  J.  Fernandez,  "Flow  Over  an  Axisvmmetric  Body  in  a  Cylindrical 
Tunnel,"  TM  79-31,  Applied  Research  Laboratory,  The  Pennsylvania 
State  University,  1  March  1979. 


14  May  1980 
CHH/pjk 


14.  V.  C.  Patel,  A.  Nakayama,  and  R.  Damian,  "Measurements  in  the  Thick 
Axisymmetric  Turbulent  Boundary  Layer  Near  the  Tail  of  a  Body  of 
Revolution,"  J.  Fluid  Mech.  J53,  345  (1974). 

15.  J.  S.  Parsons  and  R.  E.  Goodson,  "The  Optimum  Shaping  of  Axisymmetric 
Bodies  for  Minimum  Drag  in  Incompressible  Flow,"  Purdue  University 
Report  ACC-72-6,  June  1972. 

16.  V.  C.  Patel  and  Y.  T.  Lee,  "Thick  Axisymmetric  Turbulent  Boundary 
Layer  and  Near  Wake  of  a  Low-Drag  Body  of  Revolution,"  Iowa 
Institute  of  Hydraulic  Research,  low a  City,  IA,  Report  No.  210, 
December  1977. 


-36- 


14  May  1980 
GHH/pjk 


Figure  1.  Flow  Field  Solution  Domains. 
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Figure  2 


Truncated  Computational  Domain  in  the  Tail-Near  Wake  Region 
of  a  Body  of  Revolution. 
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Figure  3 
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Combination. 
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Figure  4.  Body  and  Centerline  Pressure  Distribution  for  Modified 
Spheroid  -  Potential  Flow. 
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Centerline  Pressure  Distribution  for  F-57  Body  -  Potential 
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Figure  7.  Frozen  Vorticity  and  Potential  Flow  Body  and  Centerline  Pressure 
Distribution  for  F-57  Body. 
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